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We study bosonic atoms on the p-band of a two dimensional optical square lattice in the presence of 
a confining trapping potential. Using a mean-field approach, we show how the anisotropic tunneling 
for p-band particles affects the cloud of condensed atoms by characterizing the ground state density 
and the coherence properties of the atomic states both between sites and atomic flavors. In contrast 
to the usual results based on the LDA, the atomic density can become anisotropic. This anisotropic 
effect is especially pronounced in the limit of weak atom-atom interactions and of weak lattice 
amplitudes, i.e. when the properties of the ground state are mainly driven by the kinetic energies. 
We also investigate how the trap influences known properties of the non-trapped case. In particular, 
we focus on the behavior of the anti-ferromagnetic vortex-antivortex order, which for the conflned 
system, is shown to disappear at the edges of the condensed cloud. 

PACS numbers: 03.75.Lm, 03.75.Mn 



I. INTRODUCTION 

With refined experimental techniques in trapping and 
cooling, atomic gases have become prime candidates for 
studies of mesoscopic quantum phenomena [l||. Among 
different possible experimental configurations [3] , sys- 
tems of cold atoms subjected to optical lattices consti- 
tute one of the most active topics of the current research 
in the field. In the ultracold limit, these setups may 
serve as quantum simulators which can be used to test 
actual models of condensed matter theories in a precise 
way 0- In fact, the degree of experimental control in 
optical lattice systems is so great, that by tuning the 
parameters of the lattice the atoms can be moved into 
the strongly correlated regime, therefore allowing for the 
study of a variety of phenomena which include quantum 
phase transitions 0]. Beyond experimental manipula- 
tions of the ground state, the versatility of these systems 
also makes it possible to experimentally prepare certain 
excited states. In this respect, of particular interest arc 
the states of bosons restricted to the first excited energy 
bands of the lattice, the so called p-band bosons. 

Qualitatively, the physics of p-band bosons is consid- 
erably different from the well studied systems where the 
bosons are only restricted to the lowest band (s-band 
bosons). The reason for this can be intuitively under- 
stood from the isotropic square and cubic lattices, where 
the symmetry of the lattice implies a double (square lat- 
tice) and triple (cubic lattice) degeneracy 0, Q on the 
p-band. In solid state systems such degeneracies could 
be removed via Jahn- Teller effects, but since here the 
lattice is imposed from the outside, the degeneracy is ro- 
bust. This degeneracy motivates the description of the 
atomic states in terms of orbitals related to the corre- 



sponding localized Wannicr functions, characterized by a 
node in each of the spatial directions. In the direction 
of the node, the Wannicr functions are also broader and 
since this directly influences the ease of tunneling be- 
tween sites, it directly affects the dynamical properties 
of the system. Since the properties of the tunneling of 
p-band bosons are dramatically altered from the ones on 
the s-band, a rich variety of novel quantum phases 
can appear. When interactions are taken into account, 
it has also been argued that in the limit of very strong 
atom-atom interactions, atomic population can move to 
higher energy bands, affecting thus the expected ground 
state properties of ultracold atoms in optical lattices [lol - 
lisj . The broadening of the onsite wave- functions, for ex- 
ample, was experimentally verified via microwave spec- 
troscopy [l6j . In addition, signatures of (strong) inter- 
action induced higher bands physics could also be seen 
in non-equilibrium configurations, through the mapping 
of collapse-revivals structures in the atomic density [131 
(see also Ref. [Ill)- Surprising effects are also present 
in the limit of weak interactions. In fact, it was re- 
cently observed [l^, [1^ that due to unusual dispersions, 
the physics of p-band bosons appears responsible for un- 
conventional condensation, where non-zero momentum 
states [U are occupied. We should point out, however, 
that even though experiments concerning p-band physics 
have been restricted to one dimensional, square or cubic 
lattices [l^, [2^ [22], [23|, several theoretical predictions 
have been made for other lattice configurations 24 1 . 
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In experiments, optical lattice systems are generally 
subjected to an external confining trap. Although it 
is known that even for s-band bosons, the presence of 
the trap can add important features to the physics of 
the system [25|, all the aforementioned theoretical stud- 
ies of p-band bosonic systems neglect effects originating 
from the confining trap potential. Thus, it is important 
to study how the inclusion of a trap affects the p-band 
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physics. For example, in the case of a two dimensional 
(2D) lattice it is characteristic of p-band bosons to have 
tunneling coefficients with different amplitudes in differ- 
ent directions. In the non-trapped case, this property 
of anisotropic tunneling together with the properties of 
homogeneous density distributions yields a correspond- 
ing ground state which has an anti-ferromagnetic order 
with vortcx/anti- vortex states on every second site (also 
known as the state of staggered orbital angular momen- 
tum) In trapped systems, however, the property of 
anisotropic tunneling necessarily introduces density inho- 
mogeneities which break the population balance between 
different possible atomic states (here corresponding to 
the two possible orbitals of the 2D lattice). This also 
gives rise to physics beyond the one captured by using 
the local density approximation (LDA). The fate of the 
anti-ferromagnetic order in the presence of the trap is 
then unclear. 

In this paper wc study this issue and address also other 
effects and properties which arise when p-band bosons 
are confined by an external potential. Wc mostly restrict 
the analysis to 2D, but discuss how the obtained results 
generalize to 3D. The analysis is based on the ideal gas 
theory and a mean-field approach, where we assume the 
system to be deep in the region of the superfluid phase. 
We start by presenting the theoretical framework and fol- 
low with the study of the ideal system at finite tempera- 
tures, where the critical temperature for condensation in 
a non-interacting p-band bosonic gas is calculated. Wc 
then show that for a symmetric square lattice, the zero 
temperature order parameter of the condensed ground 
state is complex also in the presence of a trap, but the 
vortex/anti- vortex structure can be lost. In particular, 
the ground state for the p-band atomic densities of the 
two flavors are shown to be different except for when the 
system is driven into the Thomas-Fermi (TF) regime in 
which case we can neglect effects stemming from the ki- 
netic tunneling energy. We complete the study with an 
analysis of the zero temperature properties of an asym- 
metric lattice. We find that due to splitting of the p-band 
degeneracy, the ground state properties may be sensitive 
to small changes in the two lattice amplitudes. 

It is important to point out that our analysis is carried 
out when influence from other bands have been omitted. 
The validity of this assumption is specially tested in the 
harmonic approximation, where two p-band atoms be- 
come degenerate with one s- and one d- atom. In fact, 
due to a 'reduced final density of states for scattering pro- 
cesses' [l^l, these decays can be significantly suppressed 
and the lifetimes of the atoms in p- orbitals become 1- 
2 order of magnitude larger than typical tunneling times. 
In addition, outside the harmonic approximation as the 
case considered throughout this paper, the actual anhar- 
monicity of the lattice breaks the (p + p ^ s + d) degen- 
eracy for almost all quasi momenta, suppressing further 
such loss processes. 



II. DERIVATION OF THE EFFECTIVE MODEL 
HAMILTONIAN 

A. Hamiltonian for p-band bosons 

In terms of the field operators ^{r'), the dynamics of 
the weakly interacting Bosc gas can be described by the 
Hamiltonian 



H 



d? 



2m 



-t-l/(r^) 



(1) 



where m corresponds to the mass of the particles, C/q to 
the strength of the interparticle interaction, and V{r') 
accounts for the effects of external potentials acting 
on the system. The field operators ^(f*) and ^'^(f^) 
annihilate and create a particle at position f" respec- 
tively, and obey the standard boson commutation rela- 



tion 



*(r"),«''^(?= 



5{r" 



In this work wc con- 



sider a trapped system in 2D with l^(r') ~ Viatt{r') + 
Vtrapi'T'), where the optical lattice potential 



Viattif) = 14 sin2 (kx') + Vy sin2 {ky') 



(2) 



has amplitudes and wave vector given, respectively, by 
Va, a € and k = 27r/A, with A being the wave 

length of the applied lasers, and where 



14rap(^) 



(3) 



describes the action of an overall slowly varying harmonic 
trap with frequency Cj. 

The common practice in the study of many-body sys- 
tems subjected to periodic potentials consists in the ex- 
pansion of the many-body Hamiltonian in terms of a suit- 
able basis, generally constructed from its corresponding 
non-interacting part. In fact, the invariance under dis- 
crete translations of the lattice implies conservation of 
quasi-momcntum and an energy spectrum having a band 
structure, which therefore immediately suggest the use 
of Bloch functions. Here, however, the presence of the 
trap breaks translational invariance and implies a finite 
size for the system, consequently destroying the symme- 
tries that rigorously justify theoretical treatment in these 
terms. On the other hand, the smoothness of the poten- 
tial implies that its charac teristic length scale fulfills the 
condition I trap = ^Jh/mCj 3> A/2, and thus we can imple- 
ment the effects of the trap in each site, by only shifting 
the onsite energies and assuming that the onsite orbitals 
remain the same in the absence of a trap. This means 
that locally the system is still effectively periodic, and 
that a satisfactory approximation can be obtained from 
the traditional framework. 

Before carrying out the expansion of the field opera- 
tors we define dimcnsionless parameters by taking the 
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recoil energy Er = h^k'^/2m as the energy scale (i.e. 
all energies are scaled by this quantity) and the in- 
verse wave vector as the typical length scale I = X/2tt, 
which produces a dimensionless trap frequency given by 
w = \i'2mu} / hk^ . In these terms, the trapping poten- 
tial becomes V(f) = (^x'^ + V^) /2, with x — kx' and 
y = ky' the dimensionless positions. From now on, we 
assume these units in all the derivations so that resulting 
equations are dimensionless. As a first step, we construct 
the bosonic operators 6,yq and 6j,q which create and anni- 
hilate, respectively, one particle delocalized in the Bloch 
state (/>i/q(r) of quasi-momentum q = (g^, qy) in the i/-th 
energy band, and use it to write 

(4) 

where the j/-sum runs over all energy bands, and the 
q-sum is over the first Brillouin zone. We also use the 
above expressions to construct the site-localized Wannier 
functions, where the operators read 

Here, Rj ~ i^iiys) — i''^jx,'^jy) labels the coordinates of 
the j'th site of the lattice (j = (jxjy), jx,jy G A/"), and 
a^j (aj,j) annihilate (create) a particle in the Wannier 
state Wi/Rj(r). For completeness, the relation between 
Wannier and Bloch functions is given by 

w.^,ir}=Y,e~^^-^^cl,,^{r). (6) 
q 

As a second step in deriving an effective model described 
by the Hamiltonian of Eq. ([ij , we choose the expansion 
of the many-body Hamiltonian in terms of (O and intro- 
duce some approximations. Our option for this picture 
is justified by the fact that while considerably simpler 
for the practical implementations, the use of Wannier 
basis together with the tight-binding approximation can 
still provide a good description as long as the lattice is 
deep enough [26|. In addition to restricting the hopping 
to nearest-neighbors (tight-binding) , we truncate the ex- 
pansion of the field operators to include only the p-bands. 

As the last step of our derivation, we clarify the used 
terminology. For a square lattice, the two p-band Wan- 
nier functions at each site j are characterized by a node 
along either the x- or y- directions. Therefore we call 
atoms with orbital wavefunctions, Wxj{r) and Wyj^f), re- 
spectively as X- and y-fiavors Q, and for completeness 



we give their explicit expressions 

(7) 

Wyj{r) = Wij^{x)w2jy{y). 

From this, the nature of the node-structure becomes 
clear. It is a direct consequence of the nodal structure 
of the Wannier functions W2j{x) and 'Wij{x). An cc-flavor 
(or equivalently p^-orbital) atom, thus, not only has a 
wavefunction with a node along the x-direction, but also 
a broader distribution along x. Accordingly, the opposite 
is true for atoms in the y-fiavor. This property directly 
affects the tunneling properties of the atoms in this sys- 
tem. 

Putting everything together, we can write down the 
resulting many-body Hamiltonian 

H = Ho + Hnn + HpD, (8) 

with the ideal part given by 

Ho = - ^^tc,p a]3i«/3j + X! X! ^trap(Rj)riaj, (9) 

where refers to the sum over nearest neighbors 

in the direction a (a,/3 = x,y) and fia^ = '^Ij'^ctj the 
atom number operator; and where the interaction terms 

a j al3,a^fi j 

(10) 

and 

account, respectively, for contribution of density-density 
and interfiavor conversion interactions. The expression 
for the interaction parameters is given by 

Ucf, ^UoJ dr\w^^{f)\'\wfsi{f)\', (12) 

and for the tunneling coefficients by 

= - / drw^^ir^ [-V^ + V{r)] (f), (13) 

where by j -f 1^ we indicate the neighboring site of j in 
the direction P, and Uq = tJol^/Er, Vp = Vp/Er are 
the dimensionless interparticle strength and lattice am- 
plitudes, respectively. From here, after substitution of 
the Wannier functions ([7]) into the above equation (|13p . 
it is straightforward to see that contributions for the tun- 
neling coefficient in the direction perpendicular to the 
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node depend uniquely from Wannier functions of the first 
band (i.e. = 1), while in the direction of the node it 
solely depends on the second band Wannier functions 
= 2). As a consequence, an x-flavor atom has larger 
probability of tunneling in the cc-direction than in the 
y-direction, while the opposite also holds for a y-flavor 
atom. We continue discussions regarding the effects of 
this anisotropic tunneling in Sec. IIIIl Also, before pro- 
ceeding with the mean-field derivations, we make a brief 
comment on the symmetries of the Hamiltonian ([5|). As 
pointed out in Ref. @, this Hamiltonian has an associ- 
ated Z2 symmetry, related to the parity of atomic flavors: 
since atom scattering processes occur in pairs, the num- 
ber of a;-flavor atoms Nx and y-flavor atoms Ny are pre- 
served modulo 2. Isotropic lattices support, in addition, 
a symmetry corresponding to swapping of atomic flavors 
X -H- J/. We will also discuss how this property implies 
a double degeneracy of the ground state for the infinite 
system. 



with the mean-field Hamiltonian 

a.P (ij)c, " j 

+X X Y (^1 + yf) 

+ J2 X ^«'3'^«J'^« + X X^ 

(17) 

and where the Hamiltonian has been normally or- 
dered prior to calculation of the coherent state expecta- 
tion value. Here the density of the flavor a is given by 
Uaj = \ipaj\'^ and normalization was imposed in the whole 
lattice as 



B. Mean-field Hamiltonian 



Except when otherwise stated, all our results follow 
from analysis of the 2D lattice. We assume the conden- 
sate confined in the transverse z— direction, and thus at 
each lattice site the system could either be purely two- 
dimensional or form condensed tubes with typically a few 
hundred of atoms [27| . In either configuration, a mean- 
field treatment is expected to give a reliable picture of 
the relevant physics [9|. 

At a mean-field level, the operators a^j are replaced by 
the complex numbers ipaj- This approximation is equiv- 
alent to assigning a coherent state at each site, j^") = 
<8>j = <8>j IV'xj,V'yj)j such that dajl"^) = ipajl"^). In 
terms of the Fock basis, the single site many-body wave- 
function reads 



= cxp 



where In) 



1^: 



(14) 



WxiTiy)-} represents the state of x-flavor 
atoms and Uy y-flavor atoms at site j. Moreover, in this 
language the onsite order parameter of site j and flavor 
a reads 'ipaj = {^\aai\^)- 

With the coherent state ansatz we can obtain the equa- 
tions of motion for the order parameter Vaj from the 
Euler-Lagrangc equations 



dL d dL 
where the Lagrangian is given by 



(15) 



AIF, 



(16) 



with N accounting for the total number of atoms. 

The Euler-Lagrange equations then correspond to a 
set of coupled Gross-Pitaevskii equations, one for each 
atomic a-flavor at each site j: 



dt 



' dt 



fS(^{x,y} 



-{Uxy + Uyxmrx; 



Pe{x,y} 
^2 

+^{^l + yl)^yi 



HUyx + Uxy)i>l-^ryy 

(19) 

Like all other parameters and variables, time t is a di- 
mensionless quantity. Note also that we take all the pa- 
rameters entering the above equations from numerically 
obtained Wannier overlap integrals according to Eqs. p2)) 
and (jl3|) . and consequently no harmonic approximation 
is imposed. This avoids some qualitatively wrong con- 
clusions which can occur with the latter assumption [11 
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III. IDEAL GAS 



uum equations 



A. Ground state properties 



\txx I 



Let us first investigate some features of the system 
in tlie non-interacting case, where the free mean-field 
Hamiltonian is given by 



a.P (ij)c 



(20) 



EET(-f+^. 



In the absence of interflavor interactions, interfiavor on- 
site coherence is not estabhshed. However, within each 
flavor it is the tunnehng which determines how the phases 
of neighboring sites are related to each other. We thus 
characterize these properties for the ground-state by min- 
imizing (|20p . To this end, the expression for the onsite 
order parameters is taken as ■i/'aj = IV'ajle*'^"-', and by 
noticing that t^xi iyy < and t 



''xy, 



tyx > we obtain 



a striped order in the phase of each flavor. More ex- 
plicitly, the phase of the x-flavor order parameter can be 
expressed as = 4'x{jx,jy) = tt x mod This 
means that neighboring sites will always keep the same 
phase in the direction perpendicular of the node, while 
in the parallel direction the phase difference will be tt. 



The discrete model (|20p can in principle be solved an- 
alytically by noticing that the Hamiltonian matrix has 
the same structure as the one of the Mathieu equation 
expanded in momentum eigenstates [l^ . The solutions is 
not very instructive as it is determined from the Fourier 
expansion of the Mathieu functions, i.e. by the transfor- 
mation matrix between quasi- and real momentum. A 
simple physical picture of the influence of the trap in the 
discrete model is instead better analyzed in the contin- 
uum limit where the analytical solutions can be given in 
closed forms. Here it is convenient to work with the order 
parameters without phase modulation. We thus impose 
the correct phase imprint responsible for rendering the 
striped order into the wavefunction ansatz. Under these 
circumstances, the phase factors can be absorbed into the 
redefinition of the tunneling coefhcient, taa ~^ —taa- In 
addition, the continuum limit consists in 'i/j^.j — >■ ipaix^ y), 
and the kinetic energy transforms as 



V'aj-Hl^ - '2'lpo.i + Ipaj-l^ > ■^1paia,/3). (21) 



With this approximation, we obtain the following contin- 



i'x{x,y), 



(22) 



■4^y{x,y), 



where x and y are dimensionless. By introducing the ef- 
fective mass rriap = |tQ,/3|~"'^/2 and parallel and transverse 
frequencies 



^± = ^\J'^\'taf}\-, a = 



(23) 



Eq. ([22]) can be written as 



'i-Q^ipxix,y) 



2mxx 



-y 



■4^a{x,y), 



2 2 

(24) 

with a similar equation for the y-flavor. We find, there- 
fore, that the continuum approximation reduces the sys- 
tem to two 2D anisotropic harmonic oscillators. It is im- 
portant to stress though, that in order to derive Eq. ()24p , 
the striped order must be correctly implemented. If the 
phase modulation is not considered before imposition of 
the continuum approximation, the resulting Hamiltonian 
is not bounded from below, and since the lattice natu- 
rally introduces a momentum cut-ofF A = tt/A at the 
edges of the Brillouin zone, it is a property not present 
in the discrete model. The initial phase imprint is thus a 
tool to circumvent this problem, where the overall effect 
of the procedure translates into inversion of the p-band 
and shifting of its minimum to the center of the Brillouin 
zone. 

In the continuum model, the anisotropy arising from 
the different tunneling elements t^x and t^y is directly 
reflected in the direction-dependence of m^^g and ujap. 
Therefore, it follows from this anisotropy that the con- 
tinuum Gaussian ground state will have different widths 
in the two directions x and y. We use this fact to define 
the anisotropy parameter 



(A.2/)2 



(25) 
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with equivalent expression to the y-flavor, and where 
= — {13)1, and (•••)„ represents the ex- 

pectation value taken with respect to ipaix, y). For sym- 
metry reasons SxSy = 1 must hold, and thus we call the 
x-flavor anisotropy parameter simply by S. This defini- 
tion is general and applies to both the discrete as well 
as for the continuum limit. It can be used to derive an 
explicit expression for the latter case 



and 



\txy\ 



1/4 



1/2 



(26) 



which as expected, predicts 5 = 1 for isotropic sys- 
tems (i.e., where both directions have the same tunnel- 
ing strengths). However, generally S ^ \^ and there- 
fore it reveals the existence of narrowing in the flavor 
density along one of the directions. This anisotropy is a 
consequence of the direction-dependence of the tunneling 
tap and is a result beyond the LDA. Note furthermore 
that when atom-atom interaction has been neglected, 
and Ny are preserved quantities and the actual ground 
state of the system will be determined from the prepa- 
ration process. For non-zero atom-atom interaction, Nx 
and Ny are no longer independently preserved due to the 
term (jlip and the interaction energy is minimized with 
Nx = Ny as will be seen in the next section. Now we 
continue with further discussions upon validity and ap- 
plicability of the continuum approximation. 



B. Ideal gas at finite temperatures 

For the ideal gas system, represented by the Hamilto- 
nian (PH]). it is rather straightforward to calculate finite 
temperature effects, either from direct numerical diago- 
nalization or using the analytical solutions obtained from 
Fourier expansions of Mathieu functions. Due to dis- 
cretization of (PH)) . the eigenstates in the harmonic trap 
are not the same as the usual eigenstates of the harmonic 
oscillator. Since implications of this for the thermody- 
namics of an ideal gas are not clear, we numerically solve 
the discrete 2D and also 3D Schrodinger equations for 
the eigenstates, and use these as a basis to study Bose- 
Einstein condensation on the p-band in the presence of a 
trap. 

In the continuum limit described by Eq. (|22p , the crit- 
ical temperature for the Bose-Einstein condensation in 
the harmonic trap is well known [2^ and given by 



T, 



{2D) 



cO 



,(215) 



in 2D and in 3D by 

rr'=-S?(iv/C(3))^/\ 



(27) 



(28) 



with C(3) ~ 1.20206, and where the trapping frequencies 
are defined as averages of the effective frequencies ([23|) 
as 



4a; 



\J \txx I \^xy I 



(29) 



2\l/3 



(30) 



For bosonic gases, the number Nt of thermal (non- 
condensed) atoms follows from 

1 



Nt^Y. 



exp(/3 {En-^i))-r 



(31) 



where /3 ~ Er/kBT is the inverse (dimcnsionless) tem- 
perature and /i is the chemical potential. Together with 
the eigenenergies En obtained by solving the Schrodinger 
equation, this can be used to compute the critical tem- 
perature for condensation in our lattice model. Notice 
however, that while below the critical temperature the 
chemical potential ^ is equal to the ground state energy, 
at higher temperatures it must be determined by fixing 
the total atom number to TV. 

We compare the predictions for the critical tempera- 
ture of the continuum and lattice models in Fig. [TJ As 
is seen, the general result in both the 2D and 3D sys- 
tems, consists in a somewhat lower critical temperature 
for very small atom numbers, but substantially larger 
critical temperature for high atom numbers. Such dif- 
ference is due to different density of states between the 
lattice and the continuum models. 

In a trap, the transition to the condensed state is typ- 
ically associated with pronounced changes in the atomic 
density distribution. A broad thermal distribution above 
the critical temperature acquires a bimodal structure as 
a density peak appears in the center corresponding to the 
macroscopic occupation of the condensate ground state. 
Also, as already discussed in the previous subsection, in 
the case of trapped p-band atoms the anisotropy is a new 
feature appearing in the density distribution. Above the 
critical temperature Tc, the density distribution has the 
same width in x- and y-directions, but below Tc the con- 
densate density distribution shares the properties of the 
ground state, which is anisotropic due to different tunnel- 
ing strengths in different directions. We give an example 
of this behavior in Fig. [2] by displaying the anisotropy pa- 
rameter ((25)) as a function of temperature for 1000 atoms. 
Furthermore, in Fig.[3]we show the density (ipnij) are the 
eigenstate wavefunctions) 



ntot 



\Mi) 



exp(/3 {En - m)) - 1 



(32) 



close to Tc and at T = 0, demonstrating the appearance 
of strong anisotropy (for single flavor) at low tempera- 
tures. 



IV. INTERACTING GAS 

A. Characterizing tlie ground state 

Until now we have not considered how interactions af- 
fect the system's ground state properties. Effects stem- 
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FIG. 1. (Color online) The critical temperature for the Bose- 
Einstein condensation as a function of atom number TV in (a) 
2D system and (b) 3D system. The dashed line shows the 
result based on approximating the discrete model with a con- 
tinuum one, and the solid line displays the numerically calcu- 
lated results of the discrete model. We used the dimensionless 
trap strength = 0.001 and \txx/txy\ = 20.1 which is our 
estimate for the ratio of tunneling strengths at Vx = Vy = 17. 



ming from the tunneling part and the corresponding 
phase ordering imposed in the minimization of the mean- 
field Hamiltonian were already discussed in Sec. IIIII We 
thus complete the characterization of the ground state of 
the system by repeating this analysis to the interacting 
part of Hmf- Since neighboring sites are not coupled by 
the interaction term, it is enough to consider the energy 
contribution within only one single site. In analogous 
procedure to the one used in the aforementioned anal- 
ysis, we substitute the expression ^paj = |V'aj|e"^°j for 
the onsite order parameter of the flavor a, and the re- 
sulting density-density and interflavor conversion parts 
of the mean-field Hamiltonian follow, respectively, as 



(33) 




FIG. 2. (Color online) The anisotropy parameter S of the 2D 
density distribution as a function of the txx-scaXeA tempera- 
ture for 1000 atoms, dimensionless trap strength w^/2 = 0.001 



and potential depth Vx = Vy = 17. 



and 



H 



(j) 

FD 



Uxy ~t~ Uyx 



nv>,jpcos(2( 



(34) 



Here, the term accounting for the density-density inter- 
actions is phase independent and gives no information 
about the on-site phase ordering. However, the interfla- 
vor conversion term will explicitly depend on the phase 
difference between the x- and y-flavor order parameters, 
and accordingly, establishes an onsite interflavor phase 
locking. In fact, when Uxy^ Uyx > the onsite energy is 
minimized with (pxi ~ 'Pyj ~ ±7r/2. 

Now combining the above argument with the results of 
Sec. IIIH we obtain both the on- and inter-site full phase 
coherence of the condensate within the lattice. To this 
end we adopt the position representation of the onsite 
order parameter 

V-j (r) = iJxiWxi (r) + i^yjWyj (r) (35) 

and apply the requirements of phase locking, which yield 



V-j (^) = I i'x} \wxs{r)±i\ tpyj I (f) , 



(36) 



where the ±-sign alternates between neighboring sites. 
Note that in the absence of a trap, flipping the sign on 
all the sites gives a new configuration with exactly the 
same energy. This characteristic degeneracy, related to 
the swapping of the flavors x y, was already pointed 
out earlier. By furthermore considering the orthonormal- 
ity property of Wannier functions, J (iru)*j(r)u'^i(r) ~ 
dapSji, we interpret the onsite order parameter as a spinor 



\^xi\ 

±i\ipyi\ 



(37) 



where the spatial dependence has been absorbed into the 
basis states Wxjif) and ■Wyj{r). In particular, the length 
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(a) 




FIG. 3. (Color online) The populations per site (for a single 
atomic flavor) of the 2D Bose gas close {kBT/txx = 1) to the 
condensation critical temperature (a), and at T = (b). In 
both examples, the number of atoms is 1000, dimensionless 
trapping strength /2 = 0.001, and potential depth Vx ~ 

Vy = 17. 



of the spinor defined in this way gives the number of 
atoms at site j, i.e. iVj = + IV'yjP- For having 

the same properties as a two-level system, the spinor on- 
site order parameter can be fully characterized by the 
Bloch vector Jj = (JxjT Jyj-, Jzj), where the components 



are 

In this picture, the length of the Bloch vector |Jj| = 
A'j corresponds to the number of atoms at site j, Jzj is 
the population imbalance between the two flavors, and 
due to the specific phase locking in ([57)) . we have Jx^ = 
0. We also point out that the Bloch vector constructed 
here corresponds to a mean-field version of the Schwinger 
angular momentum representation [soj . 

While the Bloch vector contains all the information 
about the spinor order parameter ([37]) . it does not con- 
tain the full information on the spatial dependence of 
the onsite order parameter (j36p . This can be most easily 
investigated in the harmonic approximation, where the 
Wannier functions are replaced by harmonic eigenstates. 
Using this description, we have 

V'^^M = Uxo\x±i\4^yo\y]e-'^ (39) 

with a being the effective width determined from the lat- 
tice amplitude. It is clear that for \'ipx} \ = IV'yjl the above 
onsite order parameter represents a vortex/anti- vortex 

state with quantization L zjil)^^"^ {'>^) = ±?/'j (r) where 
Lzj = —id^.. This is only true, however, in the harmonic 
approximation and when J^j = 0. Beyond the harmonic 
approximation this is not strictly true even when J^j = 0. 
Nevertheless, due to the properties of the Wannier func- 
tions, Eq. ([7]), a 7r/2 phase difference between flavors im- 
plies that the condensate density vanishes at the center 
of site j and that the condensate has a vortex like singu- 
larity in it. 

B. Properties in the symmetric lattice 

In the previous subsection we introduced the quanti- 
ties characterizing the physical state within each site. For 
the global properties we use the anisotropy parameter as 
defined in Eq. (^5)) . We numerically solve Eq. by 
employing the split-operator method [3l| , which is based 
on factorization of the time-evolution operator into spa- 
tial and momentum parts. This implies that the method 
is exact only in the limit of vanishingly small time step. 
Therefore, propagation is divided into small time steps 
and we verify the numerical accuracy by varying their 
size. In order to find the ground state we propagate an 
initial trial state in imaginary time until convergence has 
been reached. It is generally seen that convergence is 
faster if we assume an initial guess with the phase or- 
dering properties discussed in the previous section. It is 
also important to notice that a poor choice for the initial 
state may result in convergence to local, but not global. 
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energy minimum. To avoid this, we compare many differ- 
ent simulations were the initial trial state has been varied 
and the one with lowest final energy is assumed to be the 
ground state. The size of the grid is taken such that the 
atomic population is approximately zero at the edge of 
the grid, and in all simulations we consider a 2D system. 
The parameters of the Hamiltonian are calculated using 
the numerically obtained Wannicr functions, and conse- 
quently we do not impose the harmonic approximation. 

We have seen that the tunneling and the onsitc inter- 
action establish a phase locking according to Eq. (p6|) . In 
a system without the external trap and [/o 7^ 0, it follows 
that Jyj/Nj will either be -1-1 or —1, and the system pos- 
sesses a checkerboard structure, i.e. an anti-ferromagnet 
state with spins alternating between pointing in the pos- 
itive or negative y-direction. The condensate will thus 
show the staggered vortex/ anti- vortex structure. Within 
the validity of the tight-binding and single-band approx- 
imations, this result is exact. However, the strict vortex 

quantization L zjtpj'^'^H^) = i?/'] (r) is only precise in 
the harmonic approximation. In the presence of the trap, 
the inhomogencitics in the density together with the tun- 
neling anisotropy typically give rise to onsite interflavor 
population imbalance, which tends to break the anti- 
fcrromagnctic order and lower the onsite angular mo- 
mentum per particle from 1, which is expected from a 
quantized vortex with angular momentum along z. 

The ground state lattice populations lipxil"^ and iV'yjP 
for a system of = Vy = 17, cj = 0.005, and UqN = 1 
are displayed in Fig. |4] (a) and (b). It is clear how 
the anisotropy manifest itself, by rendering a conden- 
sate with spatially squeezed profile. In (c) we show the 
population imbalance Jjj. As we argued above, when- 
ever Jjj 7^ the anti-ferromagnetic order is broken, and 
from the figure it is evident that this is especially true 
in the edge of the condensate. To complement the re- 
sults, we also present the corresponding Bloch vectors in 
Fig. [5] (a). Since Jxj = 0, it is enough to show the Bloch 
vector in the spin yz-planc Jj = (0, J^j, Jzj)- By calling 
the horizontal axis the y-spin direction and the vertical 
axis the z-spin direction, we see that in the center of the 
condensate, the Jyj component dominates, while at the 
edge the Bloch vector no longer points along the horizon- 
tal direction demonstrating the breakdown of the anti- 
ferromagnetic order in these regions. Thus, at the center 
of the condensate where J^j sa 0, the anti-ferromagnetic 
ordering is still present. 

In Fig. [S] we show the ground state lattice populations 
for a more strongly interacting system with UqN — 15. 
In this case interactions and trap energies are larger than 
the kinetic energy and we approach the TF regime. We 
can see how the effects of the anisotropic density are now 
smoothed and the region of the center of the trap is en- 
larged. The latter also corresponds to the region where 
non-trapped like physics actually occurs, as confirmed in 
Figl5](b), by the presence of almost horizontal Bloch vec- 
tors. This also implies that now the ferromagnetic order 
extends over more sites in the lattice. 
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FIG. 4. (Color online) Upper two plots (a) and (b) show the x- 
and j/-flavor ground state population respectively. The lower 
plot (c) gives instead the corresponding population imbalance 
Jzj. The dimensionless system parameters are Vx ~ Vy = 17, 
(jj = 0.005, and UoN = 1. (Red color indicates an excess of 
a;-flavor atoms while blue regions have an excess of y-flavor 
atoms.) 
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FIG. 5. (Color online) The Bloch vector at the different lat- 
tice sites (the a;-component is strictly zero). The j/-spin direc- 
tion has been chosen along the horizontal axis and the z-spin 
direction along the vertical axis. The length of the vector rep- 
resents the density, while the offset from the horizontal axis 
indicates breakdown of the anti-ferromagnetic order. The lat- 
tice sites are marked by black dots. The upper plot (a) gives 
the results where interaction plays a minor role, UqN = 1, 
while in (b) UqN = 15 and interaction cannot be ignored. 
The rest of the parameters are the same as for Fig. [H 



Wc complete the study of the interacting system's 
ground state in the symmetric lattice by investigating 
the behavior of the anisotropy parameter ([25]) . Here, the 
relevant question to be understood is related to charac- 
terization of S when the system undergoes a transition 
to the TF regime. When the kinetic energy becomes sup- 
pressed, the anisotropy should vanish and hence S ^ 1. 
In the lattice there are two ways of suppressing the ki- 
netic energy, either by increasing the interparticle inter- 




FIG. 6. (Color online) Plots showing the x- (a) and j/-flavor 
(b) ground state populations, respectively, for the dimen- 
sionless system parameters Vx = Vy = 17, lu = 0.005, and 
UoN = 15. Due to the larger interaction, the squeezing effect 
is not as pronounced in this case compared to Fig. U 



action strength Uq directly by making use of Feshbach 
resonances, or by considering larger potential amplitudes. 
The predicted behavior of S is shown in Fig. [T] Note, 
that increasing UqN leads to a monotonic decrease of 
S until it asymptotically reaches 1. In the other case, 
where variation oiV^Vx=Vyis considered, S also ap- 
proaches 1 asymptotically, but now the behavior is not 
monotonic. This anomalous and surprising behavior does 
not appear in the continuum approximation. It should 
be noted that the continuum limit is evaluated in the 
ideal limit of J7o = 0, and we especially have that Scon 
is not approaching 1 as ^ — > oo. In this limit, on the 
other hand, any small Uq > will imply S = 1 since 
the kinetic term is negligible compared to the interac- 
tions. For small and moderate V, the continuum result 
is found to increase monotonously for increasing val- 
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FIG. 7. (Color online) The condensate anisotropy param- 
eter S as defined in Eq. (|25|) as a function of the interac- 
tion strength UoN and the lattice amplitude V — Vx = Vy. 
Whenever the amplitude or the interaction become large, the 
squeezing approaches one and the condensate enter into the 
TF regime. The dimensionless trap frequency oj — 0.005. 

ues of V. This behavior is not found for the discrete 
model, even for Uq = 0. Thus, for large amplitudes the 
discrete and continuum models predict qualitatively dif- 
ferent results for the squeezed profile of condensate in 
terms of the anisotropy parameter. We should also point 
out that for small amplitudes, typically V < 5Er [2g|, the 
tight-binding approximations break down and the results 
should not be taken too literally in this regime. 

C. Properties in the anisotropic lattice 

Asymmetry in the lattice breaks the degeneracy of 
X and y flavors. In order to investigate the effect of 
anisotropics we introduce the asymmetry parameter 



which controls the ratio between the lattice depths, such 
that R = 1 represents the symmetric lattice configura- 
tion we discussed earlier. We have numerically verified 
that the dominant effect of the asymmetry is to shift the 
energy levels of x- and y-flavors. By considering only a 
single site first, we note that in the harmonic approxima- 
tion this shift equals 

A = Ey- E,, = 2^(yR-iy (41) 

where E^ and Ey are the energies of the onsite flavors, 

i.e. Ea ^ J rf™*j(j) [-V + V)att(r)] Waj(j) and where 
the j dependence vanishes. In this single site picture, 
this splitting will have only a small effect if it is much 




1.0002 



R 



FIG. 8. (Color online) The parameter Jz as a function of 
the lattice asymmetry parameter R, for three different trap- 
ping frequencies, cj — 0.003 (red solid line), lj — 0.005 (black 
dashed line), and uj = 0.007 (blue dotted line). The vertical 
dashed thin lines indicate the typical sizes of S which deter- 
mines the transition region where the two atomic flavors co- 
exist. It is clear how S is decreased when the trap is "opened 
up" (decreasing w). The remaining dimensionless parameters 
are ?7o-^ = 1 and Vx = 17 (meaning that Vy — 17 /R). 



smaller than the characteristic interaction energy scale 

The picture becomes more complicated when we con- 
sider more sites. It can be, for example, that the region S 
around i? = 1 in which interaction mixes the two flavors 
changes as the trapping strength is varied, and in partic- 
ular, if S is small, the properties of the ground state may 
change dramatically with small variations in the various 
lattice parameters. On the other hand, if these parame- 
ters can be controlled, the physics around the degeneracy 
point might lead to novel physics similar to the adia- 
batic driving considered recently in Ref. js^ . However, 
it is worth pointing out that the present model possesses 
an additional property, namely that the x- and y-flavor 
densities are spatially different and adiabatic driving be- 
tween the two might therefore lead to macroscopic par- 
ticle flow within the trap. While interesting, this time- 
dependent aspect will be addressed elsewhere. 

The asymmetries for our square lattice can in princi- 
ple be implemented in two ways, either by considering a 
lattice with different wave vectors and ky or different 
amplitudes Vx and Vy . Here we characterize the behavior 
of the system in the latter process. The sensitivity to R 
can be analyzed, for example, in the value of the mean 
population inversion 

j 

If Jz = —1; the system consists of only y-flavor atoms, 
and Jz = +1 represents only atoms in the a;-flavor. Thus, 
Jz gives a measure of how much interaction mixes the two 
flavors. In the vicinity of i? = 1, the properties of Jz are 
illustrated in Fig. ([S]) . It clearly shows uniquely occupied 
flavors in both regions where R < 1 and R > 1. Also, as 



12 



expected, the exact point i? = 1 is characterized by equal 
sharing of population among the two flavors, and there- 
fore one recovers the properties of the degenerate system. 
As pointed out above, the non-zero interaction (Uq ^ 0) 
is crucial in order to stabilize the equal population at 

The confinement imposed by the harmonic trap im- 
plies that we are dealing with a finite size system. The 
frequency uj sets, in some sense, the system size and, as 
we discussed above, it is interesting to understand how 5 
depends on the system size. Figure ([8]) depicts the vari- 
ations of Jz around R = 1, and it is seen how these be- 
come more dramatic when ui is decreased. More precisely, 
there seems to be a one-to-one correspondence between 
the range 6 in which \Jz\ < 1 and oj, and as w — >■ the plot 
indicates also that J — 0. This suggests (for very weakly 
interacting systems) similar behavior to the one gener- 
ally exhibited by systems undergoing a first order phase 
transition (ssj . In addition, we also studied the ground 
state energy Eo{R) and found that dEo{R)/dR shows a 
pronounced change around 7? = 1 as w is decreased. We 
have also numerically verified that the range S grows for 
increasing interaction strength UqN in agreement with 
our earlier argument that interaction mixes the flavors. 

The above findings suggest that for weak interactions 
a careful adjustment of the lattice is required in order 
to study the anti-ferromagnetic properties. As interac- 
tions become stronger the anti-ferromagnetic properties 
become more robust. In experimental realizations even 
a small temperature might actually help to establish a 
phase coherence between x- and y-fiavor atoms since the 
energy gap between the ground and first excited ener- 
gies greatly decreases around the R — 1 point and in 
its vicinity one may expect population also of the first 
excited state. We furthermore notice that for non-zero 
w, the transition from one to the other extreme of Jz is 
smooth, and therefore by controlling the lattice ampli- 
tudes the system could be considered for studies of the 
many-body Landau-Zener transition [s^l or the Kibble- 
Zurek mechanism 13511. 



V. CONCLUSION 

We have investigated how a confining potential affects 
the properties of bosonic atoms residing on the p-bands 
of optical lattices. Our focus was on the 2D square lattice 
with equal lattice amplitudes in the two directions and 
we restricted our analysis to a mean-field approach. It is 
known that for a p-band square lattice model, even at a 
mean-field level the ground state forms non-trivial states 
in terms of an anti-ferromagnetic order As a re- 

sult of the anisotropic tunneling on the p-band together 
with the confinement introduced by the trap, we showed 
that the anti-ferromagnetic structure is destroyed in the 
edges of the condensate. The effects of the tunneling 
anisotropy are also manifest in the density profile of the 
atomic cloud, yielding a spatially elongated condensate 



in one of the two spatial directions, despite the isotropic 
trap. We showed how this narrowing is suppressed when 
the kinetic energy is lowered, either due to increasing of 
the strength of atom-atom interactions and/or due to in- 
creasing the lattice amplitudes. The same suppression 
was found also for the ideal gas when the temperature 
is increased and thereby the properties of the gas are 
greatly determined by thermal atoms. By considering 
unequal lattice amplitudes in the x- and y-directions, 
the degeneracy on the p-bands is broken, and we demon- 
strated that the sensitivity of the ground state properties 
depend strongly on the system "size". The results pre- 
sented are for 2D lattices, but it is understood that the 
general findings directly generalize to 3D as well. In the 
3D cubic case, the phase ordering can be more compli- 
cated 0, but as in the 2D case, this ordering would also 
be destroyed in the edges of the condensate in a trapped 
system. 

One point we have not addressed concerns experimen- 
tal realizations. The main source for dissipation and de- 
cohcrence in the square lattices is scattering of two p- 
band atoms into one s- and one c?-band atom [23| . 
This process is resonant in the harmonic approxima- 
tion, while it is generally off-resonant for actual lattices, 
which causes the typical life-time for p-band atoms to be 
much larger than the characteristic tunneling times. In 
Ref. [2^, coherence of p-band atoms in a cubic lattice 
was indeed demonstrated. Alternatives for suppressing 
this decay further include loading fermionic atoms into 
the s-band of the lattice [36j or consideriiig experimental 
setups with non-separable lattices [ll,[23,[33 ■ In the first 
case, the presence of fermions in the s-band prevents the 
bosonic p- band atoms to occupy the lowest band due to 
atom-atom interactions. Now in configurations involving 
non-separable lattices (e. g. superlattices), few bands 
can be separated from the rest, and thus the role of the 
{p + p — >■ s + d) scattering becomes overshadowed. In 
Refs. 0, [13], however, the experimental setup gives rise 
to hybridization of different flavor atoms and the anal- 
ysis becomes more complex than the one for the simple 
square lattice considered here. 

Another important experimentally relevant question 
concerns detection of the presented predictions. If the 
detection makes no difference between x- and y-flavor 
atoms, the Bloch vector cannot be fully measured. How- 
ever, in a recent work it was suggested how such mea- 
surements can indeed be performed [s^. The idea uti- 
lizes Raman pulses that rotate the spinor ( [57)1 similar to 
qubit measurements in atomic physics 139 11. Moreover, 



in a recent experiment on triangular lattices j27| it was 
demonstrated how the phase of the condensate affects 
the densities in time-of-flight measurements. We have 
numerically studied the full condensate order parameter 
^(x, y), and found that coherence within single sites are 
seen in ^(x,y) while long range coherence is manifested 
in the momentum distribution of ^'(a;,?/). This means 
that if the condensate density is detected at 

different time instants in a time-of-flight measurements. 
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one could in principle extract all information about the 
phase coherence. 

We believe that entering the more strongly correlated 
regime where quantum fluctuations become more im- 
portant would be of interest. The mean-field method 
adopted here is not capable of capturing these effects, 
and we therefore leave this investigation for the future. 
We esTCcially intend to study the "wedding cake" struc- 
ture [25| formed by alternating insulating Motts and su- 
perfluids in the presence of a harmonic trap, as well as 
non-equilibrium properties of the system. 
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